#ifndef POISSON_H
#define POISSON_H

// #include <cblas.h>
// #include <lapacke.h>
#include <stdexcept>
#include <cmath> 
#include <vector>
#include <iostream>
#include <memory>
#include <Eigen/Sparse>
#include "Function.hpp"

// Problem domain and boundary conditions configuration
struct DomainParams {
    enum Type { SQUARE, SQUARE_MINUS_DISK };
    Type type;
    double Point_x, Point_y, radius; // For disk
    int grid_resolution;   // M
};

struct BoundaryConditions {
    enum Type { DIRICHLET, NEUMANN, Mixed };
    Type outer_bc, inner_bc;

    std::vector<int> flag_outer_u;
    std::vector<int> one_point_loc;
    double one_point_value;
    std::vector<double> value;
    std::shared_ptr<Function> u;
    std::shared_ptr<Function> u_x;
    std::shared_ptr<Function> u_y;
    std::shared_ptr<Function> u_1;
    std::shared_ptr<Function> u_2;
    std::shared_ptr<Function> u_3;
    std::shared_ptr<Function> u_4;

    std::shared_ptr<Function> inner_u;
    std::shared_ptr<Function> inner_u_x;
    std::shared_ptr<Function> inner_u_y;
    int flag_inner_u;
};

class PoissonSolver {
    private:
        DomainParams domain;
        BoundaryConditions bc;
        double h;
        Eigen::SparseMatrix<double> A;
        std::vector<double> F;
        std::vector<double> U;
        std::shared_ptr<Function> source_term;
        std::vector<int> ipiv;
        std::vector<std::vector<int>> index_map;

        bool in_square(double x, double y) const {
            return (x >= 0 && x <= 1 && y >= 0 && y <= 1);
        }

        bool in_disk(double x, double y) const {
            double dx = x - domain.Point_x;
            double dy = y - domain.Point_y;
            return ((dx*dx + dy*dy) <= domain.radius*domain.radius);
        }

        // Check if point (x,y) is inside the domain
        bool isInside(double x, double y) const {
            if (domain.type == DomainParams::SQUARE) {
                return (x >= 0 && x <= 1 && y >= 0 && y <= 1);
            } else {
                return in_square(x,y) && !in_disk(x,y);
            }
        }

        // Verify the neighbour of the irregular point
        std::vector<double> Verify_Point(int i, int j, int category) {
            std::vector<double> q_1(2,0.0);
            std::vector<double> q_2(2,0.0);
            std::vector<double> p(2,0.0);
            p[0] = i*h;
            p[1] = j*h;

            // Horizontal direction (category == 1)
            if (category == 1) {
                q_1[1] = p[1];
                q_2[1] = p[1];
                q_1[0] = domain.Point_x - std::sqrt(domain.radius * domain.radius - (domain.Point_y - j * h) * (domain.Point_y - j * h));
                q_2[0] = domain.Point_x + std::sqrt(domain.radius * domain.radius - (domain.Point_y - j * h) * (domain.Point_y - j * h));
                if (abs(q_1[0] - p[0]) >= abs(q_2[0] - p[0])) return q_2;
                else return q_1;
            }
            
            // Vertical direction (category == 0)
            if (category == 0) {
                q_1[0] = p[0];
                q_2[0] = p[0];
                q_1[1] = domain.Point_y - std::sqrt(domain.radius * domain.radius - (domain.Point_x - i * h) * (domain.Point_x - i * h));
                q_2[1] = domain.Point_y + std::sqrt(domain.radius * domain.radius - (domain.Point_x - i * h) * (domain.Point_x - i * h));
                if (abs(q_1[1] - p[1]) >= abs(q_2[1] - p[1])) return q_2;
                else return q_1;
            }
            return q_1;
        }
        
        // Validate disk parameters
        bool validateDisk() const {
            if (domain.type != DomainParams::SQUARE_MINUS_DISK) return true;

            // Check if disk is within the square and covers ≥4 points
            double max_dist = 0.0;
            std::vector<std::pair<double, double>> corners = {{0,0}, {0,1}, {1,0}, {1,1}};
            for (const auto& corner : corners) {
                double dx = corner.first - domain.Point_x;
                double dy = corner.second - domain.Point_y;
                max_dist = std::max(max_dist, std::sqrt(dx*dx + dy*dy));
            }
            if (domain.radius >= max_dist) {
                std::cerr << "Error: Disk covers the entire square.\n";
                return false;
            }
            // Check grid points inside disk
            int count = 0;
            double h_1 = 1.0 / (domain.grid_resolution + 1);
            for (int i = 1; i <= domain.grid_resolution; ++i) {
                for (int j = 1; j <= domain.grid_resolution; ++j) {
                    double x = i * h_1;
                    double y = j * h_1;
                    if (!isInside(x, y)) count++;
                }
                if (count >= 4) break;
            }
            if (count < 4) {
                std::cerr << "Error: Disk covers <4 grid points.\n";
                return false;
            }
            return true;
        }

    public:
        PoissonSolver(const DomainParams& dp, const BoundaryConditions& bc, std::shared_ptr<Function> f)
        : domain(dp), bc(bc), source_term(f) {
            if (domain.type == DomainParams::SQUARE_MINUS_DISK && !validateDisk()) {
                std::cerr << "Error: Invaild Domain.\n";
                exit(1);
            }
        }
        
        void Build_System() {
            h = 1.0 / (domain.grid_resolution + 1);
            int m = domain.grid_resolution;
            int total_points = (m+2)*(m+2);
            int num = 0;
            double h2 = h*h;
            int flag = 0;
            int count = 0;
            index_map.resize( m+2, std::vector<int>(m+2, -1) );
            
            if (domain.type == DomainParams::SQUARE) {
                A.resize(total_points, total_points);  // 设置稀疏矩阵的大小
                A.reserve(Eigen::VectorXi::Constant(total_points, 5));  // 预分配每行最多 5 个非零元素
                U.resize(total_points, 0.0);
                F.resize(total_points, 0.0);
                ipiv.resize(total_points, 0);
            }
            
            if (domain.type == DomainParams::SQUARE && bc.outer_bc == BoundaryConditions::DIRICHLET) {
                for (int i = 0; i < total_points; ++i) {
                    flag = i % (m + 2);
                    if ( i <= m+1 ) {
                        A.insert(i,i) = 1;
                    } else if ( i >= total_points-m-2 ) {
                        A.insert(i,i) = 1;
                    } else if ( flag == 0 || flag == m + 1 ) {
                        A.insert(i,i) = 1;
                    } else { 
                        A.insert(i,i) = 4.0/h2;
                        A.insert(i,i-1) = A.insert(i,i+1) = -1.0/h2;
                        A.insert(i,i-m-2) = A.insert(i,i+m+2) = -1.0/h2;
                    }
                    
                }

                /*for (int i = 0; i < total_points; ++i) {
                    for (int j = 0; j < total_points; ++j ) {
                        std::cout<< A[i][j] <<" ";
                    }
                    std::cout<<std::endl;
                }*/
                
                for (int j = 0; j < m+2; ++j) {
                    for (int i = 0; i < m+2; ++i) {
                        if (j == 0) {
                            F[i+j*(m+2)] = bc.u_1->operator()(i*h,j*h);
                            ++ count;
                        }else if (i== m+1) {
                            F[i+j*(m+2)] = bc.u_2->operator()(i*h,j*h);
                            ++ count;
                        } else if (j == m+1) {
                            F[i+j*(m+2)] = bc.u_3->operator()(i*h,j*h);
                            ++ count;
                        } else if (i == 0) {
                            F[i+j*(m+2)] = bc.u_4->operator()(i*h,j*h);
                            ++ count;
                        } else {
                            F[i+j*(m+2)] = (*source_term)(i*h,j*h);
                        }
                    }
                }
            }
 
            if (domain.type == DomainParams::SQUARE && bc.outer_bc == BoundaryConditions::NEUMANN) {
                for (int i = 0; i < total_points; ++i) {
                    flag = i % (m + 2);
                    if (i == bc.one_point_loc[0]+(m+2)*bc.one_point_loc[1]) {
                        A.insert(i,i) = 1;
                    } else if ( i == 0 ) {
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i+1) = -1.0/h;
                        A.insert(i,i+m+2) = -1.0/h;
                    } else if ( i == m+1 ) {
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i-1) = -1.0/h;
                        A.insert(i,i+m+2) = -1.0/h;
                    } else if ( i == total_points-m-2 ) {
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i+1) = -1.0/h;
                        A.insert(i,i-m-2) = -1.0/h;
                    } else if ( i == total_points-1 ) {
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i-1) = -1.0/h;
                        A.insert(i,i-m-2) = -1.0/h;
                    } else if ( ( i >= 1 && i <= m ) ) {
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i-1) = -1.0/(2*h);
                        A.insert(i,i+1) = -1.0/(2*h);
                        A.insert(i,i+m+2) = -2.0/(2*h);
                    } else if ( ( i >= total_points-m-1 && i <= total_points -2 ) ) {
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i+1) = -1.0/(2*h);
                        A.insert(i,i-1) = -1.0/(2*h);
                        A.insert(i,i-m-2) = -2.0/(2*h);
                    } else if ( flag == 0 ) { 
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i+1) = -2.0/(2*h);
                        A.insert(i,i+m+2) = -1.0/(2*h);
                        A.insert(i,i-m-2) = -1.0/(2*h);
                    } else if ( flag == m+1 ) { 
                        A.insert(i,i) = 2.0/h;
                        A.insert(i,i-1) = -1.0/h;
                        A.insert(i,i+m+2) = -1.0/(2*h);
                        A.insert(i,i-m-2) = -1.0/(2*h);
                    } else { 
                        A.insert(i,i) = 4.0/h2;
                        A.insert(i,i+1) = A.insert(i,i-1) = -1.0/h2;
                        A.insert(i,i+m+2) = A.insert(i,i-m-2) = -1.0/h2;
                    }
                }

                for (int j = 0; j < m+2; ++j) {
                    for (int i = 0; i < m+2; ++i) {
                        if (i == bc.one_point_loc[0] && j == bc.one_point_loc[1]) {
                            F[i+j*(m+2)] = bc.one_point_value;
                        }else if (i == 0 && j == 0) {
                            F[i+j*(m+2)] = -bc.u_x->operator()(i*h,j*h) - bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        }
                        else if (j == 0 && i == m+1) {
                            F[i+j*(m+2)] = bc.u_x->operator()(i*h,j*h) - bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        }
                        else if (j == m+1 && i == 0) {
                            F[i+j*(m+2)] = -bc.u_x->operator()(i*h,j*h) + bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        }
                        else if (j == m+1 && i == m+1) {
                            F[i+j*(m+2)] = bc.u_y->operator()(i*h,j*h) + bc.u_x->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        }
                        else if (j == 0 && i >= 1 && i <= m) {
                            F[i+j*(m+2)] = - bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        }
                        else if (i == m+1 && j >= 1 && j <= m) {
                            F[i+j*(m+2)] =  bc.u_x->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        }
                        else if (j == m+1 && i >= 1 && i <= m) {
                            F[i+j*(m+2)] =  bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        }
                        else if (i == 0 && j >= 1 && j <= m) {
                            F[i+j*(m+2)] = -bc.u_x->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        } 
                        else {
                            F[i+j*(m+2)] = (*source_term)(i*h,j*h);
                        }
                    }
                }
                
            }
            
            if (domain.type == DomainParams::SQUARE && bc.outer_bc == BoundaryConditions::Mixed) {
                for (int i = 0; i < total_points; ++i) {
                    flag = i % (m + 2);
                    if ( i == 0 ) {
                            if (bc.flag_outer_u[0] == 0 && bc.flag_outer_u[3] == 0) {
                                A.insert(i,i) = 1;
                            } else if (bc.flag_outer_u[0] == 0 && bc.flag_outer_u[3] == 1) {
                                A.insert(i,i) = 1;
                            } else if (bc.flag_outer_u[0] == 1 && bc.flag_outer_u[3] == 1) {
                                A.insert(i,i) = 2.0/h;
                                A.insert(i,i+1) = -1.0/h;
                                A.insert(i,i+m+2) = -1.0/h;
                            } else if (bc.flag_outer_u[0] == 1 && bc.flag_outer_u[3] == 0) {
                                A.insert(i,i) = 1;
                            }
                    } else if ( i == m+1 ) {
                        if (bc.flag_outer_u[0] == 0 && bc.flag_outer_u[1] == 0) {
                            A.insert(i,i) = 1;
                        } else if (bc.flag_outer_u[0] == 0 && bc.flag_outer_u[1] == 1) {
                            A.insert(i,i) = 1;
                        } else if (bc.flag_outer_u[0] == 1 && bc.flag_outer_u[1] == 1) {
                            A.insert(i,i) = 2.0/h;
                            A.insert(i,i-1) = -1.0/h;
                            A.insert(i,i+m+2) = -1.0/h;
                        } else if (bc.flag_outer_u[0] == 1 && bc.flag_outer_u[1] == 0) {
                            A.insert(i,i) = 1;
                        }
                    } else if ( i == total_points-m-2 ) {
                        if (bc.flag_outer_u[2] == 0 && bc.flag_outer_u[3] == 0) {
                            A.insert(i,i) = 1;
                        } else if (bc.flag_outer_u[2] == 0 && bc.flag_outer_u[3] == 1) {
                            A.insert(i,i) = 1;
                        } else if (bc.flag_outer_u[2] == 1 && bc.flag_outer_u[3] == 1) {
                            A.insert(i,i) = 2.0/h;
                            A.insert(i,i+1) = -1.0/h;
                            A.insert(i,i-m-2) = -1.0/h;
                        } else if (bc.flag_outer_u[2] == 1 && bc.flag_outer_u[3] == 0) {
                            A.insert(i,i) = 1;
                        }
                    } else if ( i == total_points - 1 ) {
                        if (bc.flag_outer_u[1] == 0 && bc.flag_outer_u[2] == 0) {
                            A.insert(i,i) = 1;
                        } else if (bc.flag_outer_u[1] == 0 && bc.flag_outer_u[2] == 1) {
                            A.insert(i,i) = 1;
                        } else if (bc.flag_outer_u[1] == 1 && bc.flag_outer_u[2] == 1) {
                            A.insert(i,i) = 2.0/h;
                            A.insert(i,i-1) = -1.0/h;
                            A.insert(i,i-m-2) = -1.0/h;
                        } else if (bc.flag_outer_u[1] == 1 && bc.flag_outer_u[2] == 0) {
                            A.insert(i,i) = 1;
                        }
                    } else if ( i >= 1 && i <= m ) {
                        if ( bc.flag_outer_u[0] == 0 ) {
                            A.insert(i,i) = 1;
                        } else {
                            A.insert(i,i-1) = -1.0/(2*h);
                            A.insert(i,i+1) = -1.0/(2*h);
                            A.insert(i,i) = 2.0/h;
                            A.insert(i,i+m+2) = -1.0/h;
                        }
                    } else if ( i >= total_points-m-1 && i <= total_points -2 ) {
                        if ( bc.flag_outer_u[2] == 0 ) {
                            A.insert(i,i) = 1;
                        } else {
                            A.insert(i,i-1) = -1.0/(2*h);
                            A.insert(i,i+1) = -1.0/(2*h);
                            A.insert(i,i) = 2.0/h;
                            A.insert(i,i-m-2) = -1.0/h;
                        }
                    } else if ( flag == 0 && (i/(m+2)) >= 1 && (i/(m+2)) <= m ) { 
                        if ( bc.flag_outer_u[3] == 0 ) {
                            A.insert(i,i) = 1;
                        } else {
                            A.insert(i,i-m-2) = -1.0/(2*h);
                            A.insert(i,i+m+2) = -1.0/(2*h);
                            A.insert(i,i) = 2.0/h;
                            A.insert(i,i+1) = -1.0/h;
                        }
                    } else if ( flag == m+1 && (i/(m+2)) >= 1 && (i/(m+2)) <= m) { 
                        if ( bc.flag_outer_u[1]  == 0 ){
                            A.insert(i,i) = 1;
                        } else {
                            A.insert(i,i-m-2) = -1.0/(2*h);
                            A.insert(i,i+m+2) = -1.0/(2*h);
                            A.insert(i,i) = 2.0/h;
                            A.insert(i,i-1) = -1.0/h;
                        }   
                    } else { 
                        A.insert(i,i) = 4.0/h2;
                        A.insert(i,i+1) = A.insert(i,i-1) = -1.0/h2;
                        A.insert(i,i+m+2) = A.insert(i,i-m-2) = -1.0/h2;
                    }
                }
                
                for (int j = 0; j < m+2; ++j) {
                    for (int i = 0; i < m+2; ++i) {
                        if (j == 0 && i == 0) {
                            if (bc.flag_outer_u[0] == 0) {
                                F[i+j*(m+2)] = bc.u_1->operator()(i*h,j*h);
                            } else if (bc.flag_outer_u[3] == 0) {
                                F[i+j*(m+2)] = bc.u_4->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) - bc.u_x->operator()(i*h,j*h) - bc.u_y->operator()(i*h,j*h);
                            }
                        } else if (i == m+1 && j == 0) {
                            if (bc.flag_outer_u[0] == 0) {
                                F[i+j*(m+2)] = bc.u_1->operator()(i*h,j*h);
                            } else if (bc.flag_outer_u[1] == 0) {
                                F[i+j*(m+2)] = bc.u_2->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) + bc.u_x->operator()(i*h,j*h) - bc.u_y->operator()(i*h,j*h);
                            }
                        } else if (j == m+1 && i == 0) {
                            if (bc.flag_outer_u[3] == 0) {
                                F[i+j*(m+2)] = bc.u_4->operator()(i*h,j*h);
                            } else if (bc.flag_outer_u[2] == 0) {
                                F[i+j*(m+2)] = bc.u_3->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) - bc.u_x->operator()(i*h,j*h) + bc.u_y->operator()(i*h,j*h);
                            }
                        } else if (i == m+1 && j == m+1) {
                            if (bc.flag_outer_u[1] == 0) {
                                F[i+j*(m+2)] = bc.u_2->operator()(i*h,j*h);
                            } else if (bc.flag_outer_u[2] == 0) {
                                F[i+j*(m+2)] = bc.u_3->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) + bc.u_x->operator()(i*h,j*h) + bc.u_y->operator()(i*h,j*h);
                            }
                        } else if (i >= 1 && i <= m && j == 0) {
                            if (bc.flag_outer_u[0] == 0) {
                                F[i+j*(m+2)] = bc.u_1->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) - bc.u_y->operator()(i*h,j*h);
                            }
                        } else if (i >= 1 && i <= m && j == m+1) {
                            if (bc.flag_outer_u[2] == 0) {
                                F[i+j*(m+2)] = bc.u_3->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) + bc.u_y->operator()(i*h,j*h);
                            }
                        } else if (j >= 1 && j <= m && i == 0) {
                            if (bc.flag_outer_u[3] == 0) {
                                F[i+j*(m+2)] = bc.u_4->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) - bc.u_x->operator()(i*h,j*h);
                            }
                        } else if (j >= 1 && j <= m && i == m+1) {
                            if (bc.flag_outer_u[1] == 0) {
                                F[i+j*(m+2)] = bc.u_2->operator()(i*h,j*h);
                            } else {
                                F[i+j*(m+2)] = (h/2)*(*source_term)(i*h,j*h) + bc.u_x->operator()(i*h,j*h);
                            }
                        } else {
                            F[i+j*(m+2)] = (*source_term)(i*h,j*h);
                        }
                    }
                }
            }
        
            if (domain.type == DomainParams::SQUARE_MINUS_DISK) {
                for (int i = 0; i <= m+1; ++i) {
                    for (int j = 0; j <= m+1; ++j) {
                        if (in_square(i*h, j*h)){
                            if (in_disk(i*h, j*h)) {
                                index_map[i][j] = -1;
                            } else {
                                index_map[i][j] = num++;
                            }
                        }
                    }
                }

                A.resize(total_points, total_points);  // 设置稀疏矩阵的大小
                A.reserve(Eigen::VectorXi::Constant(total_points, 5));  // 预分配每行最多 5 个非零元素
                U.resize(total_points, 0.0);
                F.resize(total_points, 0.0);
                ipiv.resize(total_points, 0);
            }

            if (domain.type == DomainParams::SQUARE_MINUS_DISK && bc.outer_bc == BoundaryConditions::DIRICHLET) {
                int idx = 0;
                for (int i = 0; i < m+2; ++i) {
                    for (int j = 0; j < m+2; ++j) {
                        idx = i + j*(m+2);
                        if ( index_map[i][j] == -1 ){
                            A.insert(idx,idx) = 1.0;
                            F[idx] = -1.0;
                        } else {
                            if (j == 0) {
                                A.insert(idx,idx) = 1.0;
                                F[idx] = bc.u_1->operator()(i*h,j*h);
                            } else if (i== m+1) {
                                A.insert(idx,idx) = 1.0;
                                F[idx] = bc.u_2->operator()(i*h,j*h);
                            } else if (j == m+1) {
                                A.insert(idx,idx) = 1.0;
                                F[idx] = bc.u_3->operator()(i*h,j*h);
                            } else if (i == 0) {
                                A.insert(idx,idx) = 1.0;
                                F[idx] = bc.u_4->operator()(i*h,j*h);
                            } else {
                                if ( (index_map[i-1][j] != -1) && (index_map[i+1][j] != -1) && (index_map[i][j+1] != -1) && (index_map[i][j-1] != -1) ) {
                                    A.insert(idx,idx) = 4.0/h2;
                                    A.insert(idx,idx-1) = A.insert(idx,idx+1) = -1.0/h2;
                                    A.insert(idx,idx-m-2) = A.insert(idx,idx+m+2) = -1.0/h2;
                                    F[idx] = (*source_term)(i*h,j*h);
                                } else{
                                    // F[idx] = (*source_term)(i*h,j*h);
                                    if (index_map[i+1][j] == -1) {
                                        std::vector<double> Q = Verify_Point(i,j,1);
                                        double theta = std::abs(Q[0]/h - i);
                                        // A.insert(idx,idx+1) = -2.0/(theta*(1+theta)*h2);
                                        A.insert(idx,idx-1) = -2.0/((1+theta)*h2);
                                        if (index_map[i][j+1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(P[1]/h - j);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u-> operator()((i+1)*h, j*h)))/(theta*(1+theta)*h2) + ((bc.inner_u-> operator()(i*h, (j+1)*h))*2.0)/(alpha*(1+alpha)*h2);
                                        } else if (index_map[i][j-1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(j - P[1]/h);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u-> operator()((i+1)*h, j*h)))/(theta*(1+theta)*h2) + ((bc.inner_u-> operator()(i*h, (j-1)*h))*2.0)/(alpha*(1+alpha)*h2);
                                        
                                        } else {
                                            A.insert(idx,idx) = 2/h2 + 2/(h2*theta) ; 
                                            A.insert(idx,idx-m-2) = -1.0/h2;
                                            A.insert(idx,idx+m+2) = -1.0/h2;
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u-> operator()((i+1)*h, j*h)))/(theta*(1+theta)*h2);
                                        }
                                        
                                    } else if (index_map[i-1][j] == -1) {
                                        std::vector<double> Q = Verify_Point(i,j,1);
                                        double theta = std::abs(i - Q[0]/h);
                                        // A.insert(idx,idx-1) = -2.0/(theta*(1+theta)*h2);
                                        A.insert(idx,idx+1) = -2.0/((1+theta)*h2);
                                        if (index_map[i][j+1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(P[1]/h - j);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u-> operator()((i-1)*h, j*h)))/(theta*(1+theta)*h2) + ((bc.inner_u-> operator()(i*h, (j+1)*h))*2.0)/(alpha*(1+alpha)*h2);
                                        } else if (index_map[i][j-1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(j - P[1]/h);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u-> operator()((i-1)*h, j*h)))/(theta*(1+theta)*h2) + ((bc.inner_u-> operator()(i*h, (j-1)*h))*2.0)/(alpha*(1+alpha)*h2);
                                        
                                        } else {
                                            A.insert(idx,idx) = 2/h2 + 2/(h2*theta) ; 
                                            A.insert(idx,idx-m-2) = -1.0/h2;
                                            A.insert(idx,idx+m+2) = -1.0/h2;
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u-> operator()((i-1)*h, j*h)))/(theta*(1+theta)*h2);
                                        
                                        }
                                    } else if (index_map[i-1][j] != -1 && index_map[i+1][j] != -1){
                                        if (index_map[i][j+1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(P[1]/h - j);
                                            A.insert(idx,idx+1) = -1.0/h2;
                                            A.insert(idx,idx-1) = -1.0/h2;
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2);
                                            // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + ((bc.inner_u-> operator()(i*h, (j+1)*h))*2.0)/(alpha*(1+alpha)*h2);
                                        } else if (index_map[i][j-1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(j - P[1]/h);
                                            A.insert(idx,idx+1) = -1.0/h2;
                                            A.insert(idx,idx-1) = -1.0/h2;
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2);
                                            // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + ((bc.inner_u-> operator()(i*h, (j-1)*h))*2.0)/(alpha*(1+alpha)*h2);
                                        
                                        }
                                    }
                                } 
                                
                            }
                        }
                    }
                    
                }
            }

            if (domain.type == DomainParams::SQUARE_MINUS_DISK && bc.outer_bc == BoundaryConditions::NEUMANN) {
                int idx = 0;
                for (int j = 0; j < m+2; j++) {
                    for (int i = 0; i < m+2; i++) {
                        idx = i + j*(m+2);
                        if (index_map[i][j] == -1) {
                            A.insert(idx, idx) = 1;
                            F[idx] = -1;
                        } 
                        else if (i == bc.one_point_loc[0] && j == bc.one_point_loc[1]) {
                            A.insert(idx, idx) = 1;
                            F[idx] = bc.one_point_value;
                        } 
                        else if (i == 0 && j == 0) {
                            // A.insert(idx, idx) = 1;
                            // F[idx] = bc.u->operator()(i*h,j*h);
                            A.insert(idx,idx) = 2.0/h;
                            A.insert(idx,idx + 1) = -1.0/h;
                            A.insert(idx,idx + m + 2) = -1.0/h;
                            F[i+j*(m+2)] = -bc.u_x->operator()(i*h,j*h) - bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        } else if (i == m + 1 && j == 0) {
                            // A.insert(idx, idx) = 1;
                            // F[idx] = bc.u->operator()(i*h,j*h);
                            A.insert(idx,idx) = 2.0/h;
                            A.insert(idx,idx - 1) = -1.0/h;
                            A.insert(idx,idx + m + 2) = -1.0/h;
                            F[i+j*(m+2)] = bc.u_x->operator()(i*h,j*h) - bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        } else if (i == 0 && j == m + 1) {
                            A.insert(idx,idx) = 2.0/h;
                            A.insert(idx,idx + 1) = -1.0/h;
                            A.insert(idx,idx - m - 2) = -1.0/h;
                            F[i+j*(m+2)] = -bc.u_x->operator()(i*h,j*h) + bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        } else if (i == m + 1 && j == m + 1) {
                            A.insert(idx,idx) = 2.0/h;
                            A.insert(idx,idx - 1) = -1.0/h;
                            A.insert(idx,idx - m - 2) = -1.0/h;
                            F[i+j*(m+2)] = bc.u_x->operator()(i*h,j*h) + bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        } else if (j == 0) {
                            // A.insert(idx, idx) = 1;
                            // F[idx] = bc.u->operator()(i*h,j*h);
                            A.insert(idx,idx) = 4.0/(2*h);
                            A.insert(idx,idx + 1) = -1.0/(2*h);
                            A.insert(idx,idx - 1) = -1.0/(2*h);
                            A.insert(idx,idx + m + 2) = -2.0/(2*h);
                            F[i+j*(m+2)] =  (h/2)*(*source_term)(i*h,j*h)-bc.u_y->operator()(i*h,j*h);
                        } else if (j == m + 1) {
                            // A.insert(idx, idx) = 1;
                            // F[idx] = bc.u->operator()(i*h,j*h);
                            A.insert(idx,idx) = 4.0/(2*h);
                            A.insert(idx,idx + 1) = -1.0/(2*h);
                            A.insert(idx,idx - 1) = -1.0/(2*h);
                            A.insert(idx,idx - m - 2) = -2.0/(2*h);
                            F[i+j*(m+2)] =  bc.u_y->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        } else if (i == 0) {
                            // A.insert(idx, idx) = 1;
                            // F[idx] = bc.u->operator()(i*h,j*h);
                            A.insert(idx,idx) = 4.0/(2*h);
                            A.insert(idx,idx + 1) = -2.0/(2*h);
                            A.insert(idx,idx - m - 2) = -1.0/(2*h);
                            A.insert(idx,idx + m + 2) = -1.0/(2*h);
                            F[i+j*(m+2)] =  (h/2)*(*source_term)(i*h,j*h)-bc.u_x->operator()(i*h,j*h);
                        } else if (i == m + 1) {
                            // A.insert(idx, idx) = 1;
                            // F[idx] = bc.u->operator()(i*h,j*h);
                            A.insert(idx,idx) = 4.0/(2*h);
                            A.insert(idx,idx - 1) = -2.0/(2*h);
                            A.insert(idx,idx - m - 2) = -1.0/(2*h);
                            A.insert(idx,idx + m + 2) = -1.0/(2*h);
                            F[i+j*(m+2)] =  bc.u_x->operator()(i*h,j*h) + (h/2)*(*source_term)(i*h,j*h);
                        } else {
                            if ( (index_map[i-1][j] != -1) && (index_map[i+1][j] != -1) && (index_map[i][j+1] != -1) && (index_map[i][j-1] != -1) ) {
                                A.insert(idx,idx) = 4.0/h2;
                                A.insert(idx,idx-1) = A.insert(idx,idx+1) = -1.0/h2;
                                A.insert(idx,idx-m-2) = A.insert(idx,idx+m+2) = -1.0/h2;
                                F[idx] = (*source_term)(i*h,j*h);
                            } else{
                                // F[idx] = (*source_term)(i*h,j*h);
                                if (index_map[i+1][j] == -1) {
                                    std::vector<double> Q = Verify_Point(i,j,1);
                                    double theta = std::abs(Q[0]/h - i);
                                    // A.insert(idx,idx+1) = -2.0/(theta*(1+theta)*h2);
                                    A.insert(idx,idx-1) = -2.0/((1+theta)*h2)- 2.0/(theta*(1+theta)*h2);
                                    if (index_map[i][j+1] == -1) {
                                        std::vector<double> P = Verify_Point(i,j,0);
                                        double alpha = std::abs(P[1]/h - j);
                                        A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                        // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                        A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                        F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u_x-> operator()((i+1)*h, j*h)))/(theta*h) + ((bc.inner_u_y-> operator()(i*h, (j+1)*h))*2.0)/(alpha*h);
                                    } else if (index_map[i][j-1] == -1) {
                                        std::vector<double> P = Verify_Point(i,j,0);
                                        double alpha = std::abs(j - P[1]/h);
                                        A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                        // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                        A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                        F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u_x-> operator()((i+1)*h, j*h)))/(theta*h) - ((bc.inner_u_y-> operator()(i*h, (j-1)*h))*2.0)/(alpha*h);
                                    } else {
                                        A.insert(idx,idx) = 2/h2 + 2/(h2*theta) ; 
                                        A.insert(idx,idx-m-2) = -1.0/h2;
                                        A.insert(idx,idx+m+2) = -1.0/h2;
                                        F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u_x-> operator()((i+1)*h, j*h)))/(theta*h);
                                    }
                                        
                                } else if (index_map[i-1][j] == -1) {
                                    std::vector<double> Q = Verify_Point(i,j,1);
                                    double theta = std::abs(i - Q[0]/h);
                                    // A.insert(idx,idx-1) = -2.0/(theta*(1+theta)*h2);
                                    A.insert(idx,idx+1) = -2.0/((1+theta)*h2) - 2.0/(theta*(1+theta)*h2);
                                    if (index_map[i][j+1] == -1) {
                                        std::vector<double> P = Verify_Point(i,j,0);
                                        double alpha = std::abs(P[1]/h - j);
                                        A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                        // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                        A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                        F[idx] = (*source_term)(i*h,j*h) - (2.0*(bc.inner_u_x-> operator()((i-1)*h, j*h)))/(theta*h) + ((bc.inner_u_y-> operator()(i*h, (j+1)*h))*2.0)/(alpha*h);
                                    } else if (index_map[i][j-1] == -1) {
                                        std::vector<double> P = Verify_Point(i,j,0);
                                        double alpha = std::abs(j - P[1]/h);
                                        A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                        // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                        A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                        F[idx] = (*source_term)(i*h,j*h) - (2.0*(bc.inner_u_x-> operator()((i-1)*h, j*h)))/(theta*h) - ((bc.inner_u_y-> operator()(i*h, (j-1)*h))*2.0)/(alpha*h);
                                    } else {
                                        A.insert(idx,idx) = 2/h2 + 2/(h2*theta) ; 
                                        A.insert(idx,idx-m-2) = -1.0/h2;
                                        A.insert(idx,idx+m+2) = -1.0/h2;
                                        F[idx] = (*source_term)(i*h,j*h) - (2.0*(bc.inner_u_x-> operator()((i-1)*h, j*h)))/(theta*h);
                                        
                                    }
                                } else if (index_map[i-1][j] != -1 && index_map[i+1][j] != -1){
                                    if (index_map[i][j+1] == -1) {
                                        std::vector<double> P = Verify_Point(i,j,0);
                                        double alpha = std::abs(P[1]/h - j);
                                        A.insert(idx,idx+1) = -1.0/h2;
                                        A.insert(idx,idx-1) = -1.0/h2;
                                        A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2);
                                        // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                        A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                        F[idx] = (*source_term)(i*h,j*h) + ((bc.inner_u_y-> operator()(i*h, (j+1)*h))*2.0)/(alpha*h);
                                    } else if (index_map[i][j-1] == -1) {
                                        std::vector<double> P = Verify_Point(i,j,0);
                                        double alpha = std::abs(j - P[1]/h);
                                        A.insert(idx,idx+1) = -1.0/h2;
                                        A.insert(idx,idx-1) = -1.0/h2;
                                        A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2);
                                        // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                        A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                        F[idx] = (*source_term)(i*h,j*h) - ((bc.inner_u_y-> operator()(i*h, (j-1)*h))*2.0)/(alpha*h);
                                    
                                    }
                                }
                            }
                        }
                    }
                }
            }

            if (domain.type == DomainParams::SQUARE_MINUS_DISK && bc.outer_bc == BoundaryConditions::Mixed) {
                int idx = 0;
                for (int j = 0; j < m + 2; ++j) {
                    for (int i = 0; i < m + 2; ++i) {
                        idx = i + j * (m + 2);
                        if (index_map[i][j] == -1) {
                            // Points inside the disk
                            A.insert(idx, idx) = 1.0;
                            F[idx] = -1.0;
                        } else {
                            // Points outside the disk
                            if (i == 0 && j == 0) {
                                // Bottom-left corner
                                if (bc.flag_outer_u[0] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_1->operator()(i * h, j * h);
                                } else if (bc.flag_outer_u[3] == 0) {
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_4->operator()(i * h, j * h);
                                } else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 2.0 / h;
                                    A.insert(idx, idx + 1) = -1.0 / h;
                                    A.insert(idx, idx + m + 2) = -1.0 / h;
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) - bc.u_x->operator()(i * h, j * h) - bc.u_y->operator()(i * h, j * h);
                                }
                            } else if (i == m + 1 && j == 0) {
                                // Bottom-right corner
                                if (bc.flag_outer_u[0] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_1->operator()(i * h, j * h);
                                } else if (bc.flag_outer_u[1] == 0) {
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_2->operator()(i * h, j * h);
                                }  else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 2.0 / h;
                                    A.insert(idx, idx - 1) = -1.0 / h;
                                    A.insert(idx, idx + m + 2) = -1.0 / h;
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) + bc.u_x->operator()(i * h, j * h) - bc.u_y->operator()(i * h, j * h);
                                }
                            } else if (i == 0 && j == m + 1) {
                                // Top-left corner
                                if (bc.flag_outer_u[3] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_4->operator()(i * h, j * h);
                                } else if (bc.flag_outer_u[2] == 0) {
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_3->operator()(i * h, j * h);
                                }  else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 2.0 / h;
                                    A.insert(idx, idx + 1) = -1.0 / h;
                                    A.insert(idx, idx - m - 2) = -1.0 / h;
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) - bc.u_x->operator()(i * h, j * h) + bc.u_y->operator()(i * h, j * h);
                                }
                            } else if (i == m + 1 && j == m + 1) {
                                // Top-right corner
                                if (bc.flag_outer_u[1] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_2->operator()(i * h, j * h);
                                } else if (bc.flag_outer_u[2] == 0) {
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_3->operator()(i * h, j * h);
                                }  else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 2.0 / h;
                                    A.insert(idx, idx - 1) = -1.0 / h;
                                    A.insert(idx, idx - m - 2) = -1.0 / h;
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) + bc.u_x->operator()(i * h, j * h) + bc.u_y->operator()(i * h, j * h);
                                }
                            } else if (j == 0) {
                                // Bottom boundary
                                if (bc.flag_outer_u[0] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_1->operator()(i * h, j * h);
                                } else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 4.0 / (2 * h);
                                    A.insert(idx, idx + 1) = -1.0 / (2 * h);
                                    A.insert(idx, idx - 1) = -1.0 / (2 * h);
                                    A.insert(idx, idx + m + 2) = -2.0 / (2 * h);
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) - bc.u_y->operator()(i * h, j * h);
                                }
                            } else if (j == m + 1) {
                                // Top boundary
                                if (bc.flag_outer_u[2] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_3->operator()(i * h, j * h);
                                } else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 4.0 / (2 * h);
                                    A.insert(idx, idx + 1) = -1.0 / (2 * h);
                                    A.insert(idx, idx - 1) = -1.0 / (2 * h);
                                    A.insert(idx, idx - m - 2) = -2.0 / (2 * h);
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) + bc.u_y->operator()(i * h, j * h);
                                }
                            } else if (i == 0) {
                                // Left boundary
                                if (bc.flag_outer_u[3] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_4->operator()(i * h, j * h);
                                } else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 4.0 / (2 * h);
                                    A.insert(idx, idx + 1) = -2.0 / (2 * h);
                                    A.insert(idx, idx - m - 2) = -1.0 / (2 * h);
                                    A.insert(idx, idx + m + 2) = -1.0 / (2 * h);
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) - bc.u_x->operator()(i * h, j * h);
                                }
                            } else if (i == m + 1) {
                                // Right boundary
                                if (bc.flag_outer_u[1] == 0) {
                                    // Dirichlet condition
                                    A.insert(idx, idx) = 1.0;
                                    F[idx] = bc.u_2->operator()(i * h, j * h);
                                } else {
                                    // Neumann condition
                                    A.insert(idx, idx) = 4.0 / (2 * h);
                                    A.insert(idx, idx - 1) = -2.0 / (2 * h);
                                    A.insert(idx, idx - m - 2) = -1.0 / (2 * h);
                                    A.insert(idx, idx + m + 2) = -1.0 / (2 * h);
                                    F[idx] = (h / 2) * (*source_term)(i * h, j * h) + bc.u_x->operator()(i * h, j * h);
                                }
                            } else {
                                // Interior points
                                if ((index_map[i - 1][j] != -1) && (index_map[i + 1][j] != -1) && (index_map[i][j + 1] != -1) && (index_map[i][j - 1] != -1)) {
                                    A.insert(idx, idx) = 4.0 / (h * h);
                                    A.insert(idx, idx - 1) = A.insert(idx, idx + 1) = -1.0 / (h * h);
                                    A.insert(idx, idx - m - 2) = A.insert(idx, idx + m + 2) = -1.0 / (h * h);
                                    F[idx] = (*source_term)(i * h, j * h);
                                } else if (bc.flag_inner_u == 0) {
                                    // Handle irregular points near the disk boundary
                                    if (index_map[i + 1][j] == -1) {
                                        std::vector<double> Q = Verify_Point(i, j, 1);
                                        double theta = std::abs(Q[0] / h - i);
                                        A.insert(idx, idx - 1) = -2.0 / ((1 + theta) * h * h);
                                        if (index_map[i][j + 1] == -1) {
                                            std::vector<double> P = Verify_Point(i, j, 0);
                                            double alpha = std::abs(P[1] / h - j);
                                            A.insert(idx, idx) = 2 / (h * h * alpha) + 2 / (h * h * theta);
                                            A.insert(idx, idx - m - 2) = -2.0 / ((1 + alpha) * h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + (2.0 * (bc.inner_u->operator()((i + 1) * h, j * h))) / (theta * (1 + theta) * h * h) + ((bc.inner_u->operator()(i * h, (j + 1) * h)) * 2.0) / (alpha * (1 + alpha) * h * h);
                                        } else if (index_map[i][j - 1] == -1) {
                                            std::vector<double> P = Verify_Point(i, j, 0);
                                            double alpha = std::abs(j - P[1] / h);
                                            A.insert(idx, idx) = 2 / (h * h * alpha) + 2 / (h * h * theta);
                                            A.insert(idx, idx + m + 2) = -2.0 / ((1 + alpha) * h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + (2.0 * (bc.inner_u->operator()((i + 1) * h, j * h))) / (theta * (1 + theta) * h * h) + ((bc.inner_u->operator()(i * h, (j - 1) * h)) * 2.0) / (alpha * (1 + alpha) * h * h);
                                        } else {
                                            A.insert(idx, idx) = 2 / (h * h) + 2 / (h * h * theta);
                                            A.insert(idx, idx - m - 2) = -1.0 / (h * h);
                                            A.insert(idx, idx + m + 2) = -1.0 / (h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + (2.0 * (bc.inner_u->operator()((i + 1) * h, j * h))) / (theta * (1 + theta) * h * h);
                                        }
                                    } else if (index_map[i - 1][j] == -1) {
                                        std::vector<double> Q = Verify_Point(i, j, 1);
                                        double theta = std::abs(i - Q[0] / h);
                                        A.insert(idx, idx + 1) = -2.0 / ((1 + theta) * h * h);
                                        if (index_map[i][j + 1] == -1) {
                                            std::vector<double> P = Verify_Point(i, j, 0);
                                            double alpha = std::abs(P[1] / h - j);
                                            A.insert(idx, idx) = 2 / (h * h * alpha) + 2 / (h * h * theta);
                                            A.insert(idx, idx - m - 2) = -2.0 / ((1 + alpha) * h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + (2.0 * (bc.inner_u->operator()((i - 1) * h, j * h))) / (theta * (1 + theta) * h * h) + ((bc.inner_u->operator()(i * h, (j + 1) * h)) * 2.0) / (alpha * (1 + alpha) * h * h);
                                        } else if (index_map[i][j - 1] == -1) {
                                            std::vector<double> P = Verify_Point(i, j, 0);
                                            double alpha = std::abs(j - P[1] / h);
                                            A.insert(idx, idx) = 2 / (h * h * alpha) + 2 / (h * h * theta);
                                            A.insert(idx, idx + m + 2) = -2.0 / ((1 + alpha) * h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + (2.0 * (bc.inner_u->operator()((i - 1) * h, j * h))) / (theta * (1 + theta) * h * h) + ((bc.inner_u->operator()(i * h, (j - 1) * h)) * 2.0) / (alpha * (1 + alpha) * h * h);
                                        } else {
                                            A.insert(idx, idx) = 2 / (h * h) + 2 / (h * h * theta);
                                            A.insert(idx, idx - m - 2) = -1.0 / (h * h);
                                            A.insert(idx, idx + m + 2) = -1.0 / (h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + (2.0 * (bc.inner_u->operator()((i - 1) * h, j * h))) / (theta * (1 + theta) * h * h);
                                        }
                                    } else if (index_map[i - 1][j] != -1 && index_map[i + 1][j] != -1) {
                                        if (index_map[i][j + 1] == -1) {
                                            std::vector<double> P = Verify_Point(i, j, 0);
                                            double alpha = std::abs(P[1] / h - j);
                                            A.insert(idx, idx + 1) = -1.0 / (h * h);
                                            A.insert(idx, idx - 1) = -1.0 / (h * h);
                                            A.insert(idx, idx) = 2 / (h * h * alpha) + 2 / (h * h);
                                            A.insert(idx, idx - m - 2) = -2.0 / ((1 + alpha) * h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + ((bc.inner_u->operator()(i * h, (j + 1) * h)) * 2.0) / (alpha * (1 + alpha) * h * h);
                                        } else if (index_map[i][j - 1] == -1) {
                                            std::vector<double> P = Verify_Point(i, j, 0);
                                            double alpha = std::abs(j - P[1] / h);
                                            A.insert(idx, idx + 1) = -1.0 / (h * h);
                                            A.insert(idx, idx - 1) = -1.0 / (h * h);
                                            A.insert(idx, idx) = 2 / (h * h * alpha) + 2 / (h * h);
                                            A.insert(idx, idx + m + 2) = -2.0 / ((1 + alpha) * h * h);
                                            F[idx] = (*source_term)(i * h, j * h) + ((bc.inner_u->operator()(i * h, (j - 1) * h)) * 2.0) / (alpha * (1 + alpha) * h * h);
                                        }
                                    }
                                } else if (bc.flag_inner_u == 1){
                                    // F[idx] = (*source_term)(i*h,j*h);
                                    if (index_map[i+1][j] == -1) {
                                        std::vector<double> Q = Verify_Point(i,j,1);
                                        double theta = std::abs(Q[0]/h - i);
                                        // A.insert(idx,idx+1) = -2.0/(theta*(1+theta)*h2);
                                        A.insert(idx,idx-1) = -2.0/((1+theta)*h2)- 2.0/(theta*(1+theta)*h2);
                                        if (index_map[i][j+1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(P[1]/h - j);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u_x-> operator()((i+1)*h, j*h)))/(theta*h) + ((bc.inner_u_y-> operator()(i*h, (j+1)*h))*2.0)/(alpha*h);
                                        } else if (index_map[i][j-1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(j - P[1]/h);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u_x-> operator()((i+1)*h, j*h)))/(theta*h) - ((bc.inner_u_y-> operator()(i*h, (j-1)*h))*2.0)/(alpha*h);
                                        } else {
                                            A.insert(idx,idx) = 2/h2 + 2/(h2*theta) ; 
                                            A.insert(idx,idx-m-2) = -1.0/h2;
                                            A.insert(idx,idx+m+2) = -1.0/h2;
                                            F[idx] = (*source_term)(i*h,j*h) + (2.0*(bc.inner_u_x-> operator()((i+1)*h, j*h)))/(theta*h);
                                        }
                                            
                                    } else if (index_map[i-1][j] == -1) {
                                        std::vector<double> Q = Verify_Point(i,j,1);
                                        double theta = std::abs(i - Q[0]/h);
                                        // A.insert(idx,idx-1) = -2.0/(theta*(1+theta)*h2);
                                        A.insert(idx,idx+1) = -2.0/((1+theta)*h2) - 2.0/(theta*(1+theta)*h2);
                                        if (index_map[i][j+1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(P[1]/h - j);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) - (2.0*(bc.inner_u_x-> operator()((i-1)*h, j*h)))/(theta*h) + ((bc.inner_u_y-> operator()(i*h, (j+1)*h))*2.0)/(alpha*h);
                                        } else if (index_map[i][j-1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(j - P[1]/h);
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2*theta);
                                            // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) - (2.0*(bc.inner_u_x-> operator()((i-1)*h, j*h)))/(theta*h) - ((bc.inner_u_y-> operator()(i*h, (j-1)*h))*2.0)/(alpha*h);
                                        } else {
                                            A.insert(idx,idx) = 2/h2 + 2/(h2*theta) ; 
                                            A.insert(idx,idx-m-2) = -1.0/h2;
                                            A.insert(idx,idx+m+2) = -1.0/h2;
                                            F[idx] = (*source_term)(i*h,j*h) - (2.0*(bc.inner_u_x-> operator()((i-1)*h, j*h)))/(theta*h);
                                            
                                        }
                                    } else if (index_map[i-1][j] != -1 && index_map[i+1][j] != -1){
                                        if (index_map[i][j+1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(P[1]/h - j);
                                            A.insert(idx,idx+1) = -1.0/h2;
                                            A.insert(idx,idx-1) = -1.0/h2;
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2);
                                            // A.insert(idx,idx+m+2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx-m-2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) + ((bc.inner_u_y-> operator()(i*h, (j+1)*h))*2.0)/(alpha*h);
                                        } else if (index_map[i][j-1] == -1) {
                                            std::vector<double> P = Verify_Point(i,j,0);
                                            double alpha = std::abs(j - P[1]/h);
                                            A.insert(idx,idx+1) = -1.0/h2;
                                            A.insert(idx,idx-1) = -1.0/h2;
                                            A.insert(idx,idx) = 2/(h2*alpha) + 2/(h2);
                                            // A.insert(idx,idx-m-2) = -2.0/(alpha*(1+alpha)*h2);
                                            A.insert(idx,idx+m+2) = -2.0/((1+alpha)*h2) - 2.0/(alpha*(1+alpha)*h2);
                                            F[idx] = (*source_term)(i*h,j*h) - ((bc.inner_u_y-> operator()(i*h, (j-1)*h))*2.0)/(alpha*h);
                                        
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        
        }
        
        void Solve() {
            // 使用 Eigen 的稀疏矩阵求解器
            Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
            solver.compute(A);
            if (solver.info() != Eigen::Success) {
                std::cerr << "Error: Sparse matrix decomposition failed.\n";
                exit(1);
            }
    
            Eigen::VectorXd F_eigen = Eigen::Map<Eigen::VectorXd>(F.data(), F.size());
            Eigen::VectorXd U_eigen = solver.solve(F_eigen);
            if (solver.info() != Eigen::Success) {
                std::cerr << "Error: Sparse matrix solver failed.\n";
                exit(1);
            }
    
            U.assign(U_eigen.data(), U_eigen.data() + U_eigen.size());
        }
        
        const std::vector<std::vector<int>>& GetMap() const {
            return index_map;
        }

        const std::vector<double>& GetSolution() const {
            return U;
        }
};
#endif